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The 'simulation' of sclf-gravitating isolated systems is one of the most challenging and hopefully rewarding tasks in 
relativity today. The challenge is to develop a reliable, accurate and maintainable code. The possible results obtained 
with this code may have consequences for different areas of classical relativity such as the mathematical question of 
the cosmic censorship hypothesis on the one hand and the prediction of wave forms for gravitational wave detectors on 



the other hand. There exist various approaches towards the development of such a code (see 39 for a recent review). 
Their majority is based on (variations of) the standard ADM equations [Q. Some of these numerical investigations 
have yielded quite impressive results p, |25fl. 



o 

■ One of the most successful approaches so far towards the simulation of isolated systems is based on the characteristic 
formulation of Einstein's equation as developed by Bondi B , Sachs Q , Newman and Penrose Q . This approach is 
presented by R. Bartnik and A. Norton |J, J. A. Font [|l5) and L. Lehner in this volume. 

The objective of this presentation is to discuss some aspects of the conformal approach to numerical relativity which 
is based on Friedrich's conformal field equations. The general conceptional and analytical background has been given 
' by Friedrich in this volume |2lj while the numerical efforts have been summarized in ^0|, see also Husa |$6|. Thus, 
we can take the opportunity to restrict ourselves to a more detailed discussion of some of the numerical problems that 
we encounter when we try to solve the conformal field equations. 

To a large extent the issues involved in the numerical treatment of the hyperboloidal initial value problem for the 
conformal field equations are the same as for any system of geometric PDEs which split into evolution and constraint 
equations. However, there are some crucial differences which are due to the fact that null infinity can be located on 
• the grid which we want to discuss here. 

Any numerical implementation of such a system of PDEs will have to face the following broad issues: (i) the 
construction of initial data, (ii) the design of a stable evolution scheme, (Hi) the implementation of (meaningful) 
boundary conditions, (iv) the control of the constraint violation during the evolution and, finally, (v) the search for 
'good' gauges. These subjects are by no means independent. In fact, the issue (v) is probably the most important 
O" 1 one because it influences all the others. But it is also the most obscure one because the 'goodness' of a gauge can in 
5^ | general not be inferred without looking at actual evolutions. Hence the 'invention' of good gauges is more like an art 
than a science. 

The design of a stable evolution scheme and the design of boundary conditions are to a large extent uncritical. 
The choice of boundary conditions is not as crucial for the conformal approach as it is for others. This is due to the 
fact that the boundary of the numerical grid is located in the unphysical region so that it is causally disconnected 
from the physical part of the grid. This means that the influence of the boundary cannot reach the interior of the 
space-time. Clearly, these statements refer only to the analytical time evolution while the discrete numerical evolution 
scheme will also propagate parasitic modes at higher speeds than the speed of light so that numerically there will be 
an influence of the boundary on the physical space-time. However, if we use a well-posed numerical formulation of 
the initial-boundary value problem then this influence should die out with the order of the scheme when we increase 
the accuracy. 

This implies that we have 'only' to make sure that our boundary conditions are compatible with the numerical 
evolution scheme in the sense that the resulting scheme is well-posed. There is a large literature in the numerical 



analysis community on the various issues of boundary conditions (see e.g. [14, |2S|, |47|) which can be consulted for- 
tius purpose. The method of choice for the time evolution is a 'method of lines' with fourth-order discretization in 
space and a high order Runge-Kutta method to solve the ensuing coupled system of ODEs. The problem which arises 
in connection with the boundary conditions is that of complexity. A high-order numerical evolution scheme needs 
boundary conditions discretized at the same order for the overall scheme to stay at the intended order. This means 
that one has to take time derivatives at the boundary and express them using the evolution equations in terms of 
the spatial derivatives. This is a task which can get very complicated for such complex evolution equations as those 
obtained from the conformal field equations. 

Another issue in relation to the boundary conditions is that they should be compatible with the constraints. This 
leads to the necessity to analyse the the mathematical initial-boundary value problem for the Einstein and conformal 
field equations. There is a recent mathematical result [E4| which treats the standard Einstein equations but the 
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extension to the conformal field equations is still lacking. There is also some related work on the initial-boundary 
value problem for linearized gravity (4^, Q and an implementation of constraint-preserving boundary conditions in 
spherical symmetry JTof . 

The problem of keeping the constraints satisfied during the numerical evolution is one which seems to plague most 
if not all the codes in numerical relativity. The issue here is that the numerical evolution tends to exponentially drive 
away the system from the constraint hypersurface although the analytical time evolution guarantees that the evolution 
is confined to the constraint surface for all times p7| . The question arises as to whether this exponential violation 
of the constraints is a purely numerical effect or whether it has a cause which can be understood analytically. This 
cause would have to lie in the way small deviations from zero in the constraints are propagated by the time evolution. 
Analytically this is regulated by the 'subsidiary system' of evolution equations for the constraints. Therefore, in order 
to shed some light onto this issue one should analyse the properties of that system in some detail. 

There remains the question of finding initial data for the time evolution. The problems which arise in this area are 
specific to the conformal approach and hence we take the opportunity to expose them here in more detail. 



II. THE ANDERSSON— CHRUSCIEL—FRIEDRICH PROCEDURE 



The first step in any numerical simulation is to provide initial data for the evolution equations. In the confor- 
mal approach this means that we have to determine the fields on the initial hypersurface in such a way that they 
obey the conformal constraints which have been obtained from the conformal field equations by a standard 3+1- 
decomposition |^0| . The structure of these equations is rather complicated and to this day they have resisted attempts 
to cast them into a form suitable for mathematical analysis or even to devise a reasonable numerical scheme for their 
solution. Recently, however, Butscher || has made some progress towards a reformulation of these equations. Ulti- 
mately, we should strive to obtain initial data directly from solving these conformal constraints. 

So far we obtain initial data by an implementation of a procedure originally due to Andersson, Chrusciel and 
Friedrich H and expanded on in The mathematical background and various results are described in detail by 
L. Andersson H in his contribution to the present volume. We will briefly describe the essential idea. The starting 
point for obtaining initial data in unphysical space-time (M, g a b) with this procedure is a hyperboloidal hypersurface 
E with induced metric h a b and extrinsic curvature k a b in physical space-time (M,g a b). 

For simplicity we assume that the extrinsic curvature of E is pure trace, k a b — \h a bk 1 the more general case is 
treated in Q . The vacuum momentum constraints imply that k is constant and the vacuum Hamiltonian constraint 
further implies that the scalar curvature R is also constant. Rescaling the induced metric by a constant factor we can 
assume that 

R = -6 (1) 

which in turn implies k — 3 (recall that a hyperboloidal hypersurface is characterized by the fact that it has an 
asymptotically constant negative curvature). 

The next step is to assume that the physical space-time admits a conformal extension so that we may regard 
(M,g a b) as being embedded in the unphysical space-time (M,g a b) where it admits a boundary . The two metrics 
are conformally equivalent 

g a b = ^ 2 g a b (2) 

with a conformal factor f2 which is positive on E and which vanishes on J^. The hyperboloidal hypersurface E extends 
as a space-like hypersurface E beyond J? intersecting it transversally. 

The embedding of E in M induces a (negative definite) metric h a b and extrinsic curvature k a b which are related on 
E to the physical data by the following relations 

h a b = Q 2 h a b, k a b = 7T (k a b — Q^abj ■ (3) 

Here, a = t a X7 a n is the normal derivative of SI with respect to the (g)-normal vector t a of E. This implies that the 
trace-free parts transform homogeneously so that if k a b is pure trace then also k a b is pure trace, while for the traces 
we obtain k = kQ — 3c 

Furthermore, the Ricci tensors of h a b and h a b are related by the usual conformal transformation formula (d a is the 
Levi-Civita connection of h a b) 

- d a d b n ao d c nd c n 

Rab Rab - Kb— + Zhab-^-^ (4) 
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from which we obtain the transformation of the scalar curvature 

R = Q 2 R - 4QAQ, + 6 d c fld c il. (5) 
Thus, we can express the physical Hamiltonian constraint in terms of unphysical quantities as 

fl 2 R - iflAfl + 6 d c fld c tt = -6. (6) 

Recall that the only condition, apart from smoothness, imposed upon the conformal factor is its vanishing on ^ . 
Thus, we may write f2 = oj(j) a for some number a where <fi is a strictly positive function on E so that u) vanishes on J 1 
to the same order as does fi. In this way we 'decouple' the two essential properties of f2: the localization of J? on E, 
indicated by the vanishing of u> which we consider to be given once and for all and secondly its role in relating the 
physical and unphysical metrics which is taken over by <p. 

Now it is easy to see that upon insertion of this expression for f2 into (||) one obtains an equation which contains 
up to second order derivatives of <f> and, in particular, terms which are proportional to d c 4>d c 4>. Choosing a = —2 one 
can eliminate these quadratic terms and arrive at the equation 

8w 2 A0 - %u:d c ujd c (t> + [u?R - AujAuj + Qd c uod c uj] <j) = -6</> 5 . (7) 

Since u> is fixed we may regard this equation as determining 4> and hence fi. This equation could be called the conformal 
Lichnerowicz equation because it arises from the Hamiltonian constraint in a similar way as the Lichnerowicz equation 
except that it is formulated in the unphysical space-time. On the other hand, the equation determines a conformal 
factor which fixes a specific member from a given conformal class of metrics which has constant scalar curvature. In 
the mathematical literature this problem is known as the Yamabe problem. So we call (|7]) the Lichnerowicz-Yamabe- 
Equation (LYE). We will discuss this equation in more detail later on so that we point out only its most relevant 
property now: the equation is singular on J 1 in the sense that for u> — all the terms containing derivatives of <fi vanish 
(we assume here that (f> and its derivatives remain bounded). There remains only an algebraic relation between the 
values of <f> and some geometric quantities. The reason for this non-regularity is easy to understand. If the equation 
was regular on J? then one would presumably be able to provide boundary data for </>. This, however, contradicts 
the physical picture of ^ being a universal construction in the sense that for any asymptotically flat space-time one 
obtains the 'same' J* . The degenerate character of the equation on J? leads to the question whether and under what 
conditions there exist unique solutions. This has been treated in detail in |3j] and we will return to this issue in 



Sect. Ill 



In order to find initial data for the conformal field equations we may thus proceed as follows. We start out with some 
3-dimensional manifold E carrying a (negative definite) Riemannian metric h a \>. Since the equation is conformally 
invariant it is really only the conformal class of the metric h a b which is relevant here. We prescribe the boundary 
function u> whose sole purpose is to fix the topology of the space-time we are trying to simulate. The only condition 
to be satisfied by this function is that its zero-set be an embedded 2-dimensional submanifold of E (it need not be 
connected nor have spherical topology). Next we try to find a non- vanishing solution <j) of (Q) with the given data. 
This provides a conformal factor Q = u><f)~ 2 which conformally rescales the metric h a b to a metric h a b on E (where 
O > 0) with constant negative curvature R = —6. Imposing k a b oc h a b implies k a b oc h a b and so k — const. In fact, 
choosing for k any smooth function on E and defining a = ^ktt — 1 guarantees that k = 3 so that the physical vacuum 

constraints will be satisfied on E. 

Once this conformal factor is determined we can continue to generate the other initial data, such as the Ricci- and 
Weyl tensors. In the present case with k a b being pure trace the relevant fields are the trace-free part of the Ricci 
tensor and the electric part of the Weyl tensor. They are given by 

$ab - ~ \d a dbSl ~ ^afeAfiJ , (8) 

E ab = ~ (*ab - Rab + ^h ab Rj . (9) 

The crucial point here is that the computation of these tensors involves division by £1. This raises two questions: 
under what conditions are these expressions well defined and how can we compute them accurately numerically. The 



first question has been answered in M while the second is the subject of Sect. IV 



III. THE LICHNEROWICZ— YAM ABE-EQUATION 

The degeneracy on J> of the Lichnerowicz- Yamabe-Equation ([?]) poses special problems for its numerical solution. 
This is due not so much to the solution process itself which can be handled quite nicely with the usual methods. The 
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problems arise from the special circumstances prevailing in the numerical setup. The results in jjl || apply to the 
situation where £ is a 3-dimensional manifold with boundary given by uj = 0. In this situation there will be a unique 
bounded solution of (Q) on E which extends smoothly to the boundary provided the data, i.e. the specified conformal 
metric, is smooth up to the boundary and satisfies a certain condition at the boundary whose physical meaning is the 
vanishing of the shear on <# . 

In numerical applications this situation although desirable cannot usually be achieved. This has several reasons. 
In 3D simulations one normally uses Cartesian coordinates while for asymptotically flat space-times J? has spherical 
'cuts'. This means that the zero-set of uj cannot be aligned with the grid and in particular it cannot be the boundary 
of the grid. Therefore, there will be nodes of the grid which lie outside of the physical region, beyond J* . Even if one 
decided to adapt the coordinates to J* so that uj = lies on a grid plane then this could be done only for at most 
two asymptotic regions like in the Kruskal space-time 1 . For a two black-hole space-time which has three asymptotic 
regions one of these would lie inside the grid and there would again be nodes of the grid which lie beyond J? in the 
unphysical region. Finally, even in the case of one J? and spherical coordinates one needs for numerical reasons a 
number of grid points beyond in order to implement boundary conditions for the time evolution. Thus, we have to 
face the situation where the grid not only covers the physical space-time but extends beyond <# into the unphysical 
region. 

Therefore, we have to ask what to do at the regions outside the physical space-time. There are essentially two ways 
to deal with this situation. We can either solve the equation only in the physical region, determine all the initial data 
there and then extend these data as smoothly as possible onto the unphysical regions. This is the path followed by 



Hiibner [B2L |33|, |3J, 35 . One major drawback of this possibility is that in general the conformal constraints will not 
be satisfied in the unphysical region because the extension does not respect them. Not only is this an esthetically 
unpleasing feature - after all, the unphysical region can be considered as an asymptotically flat space-time in its 
own right - but almost certainly it will also be counter-productive. Some features like e.g. radiation extraction and 
'scri- freezing' [ p"7| p"8[ rely on the constraints being satisfied in a neighbourhood of J? . Furthermore, it is quite likely 
that the fact that immediately beyond J? the constraints fail could allow constraint violating modes to penetrate into 
the physical space-time triggering the above mentioned exponential growth of the constraints there (for a discussion 
of this phenomenon in the context of the conformal field equations see |3^]). Of course, ideally J 1 should act as a 
one-way membrane for the propagation permitting outward flow but inhibiting inward flow. However, this is true for 
the (exact) continuous time evolution while for the (approximate) discrete evolution this depends very much on the 
method used to evolve from one time level to the next. It is true that the influence of the outside region upon the 
inside should die away with the order of the evolution scheme employed but ultimately only numerical tests will settle 
this point. 

So in order to avoid constraint violation in the unphysical region we are forced to solve the Lichnerowicz-Yamabe- 
Equation (LYE) on the entire numerical grid where ^ is not necessarily aligned with the grid boundary. However, 
this is not as straightforward as it may seem. As alluded to already above, the region outside J* can be regarded as 
an asymptotically flat space-time in its own right. The two space-times do not really have any relationship except for 
the superficial property that they match on the same numerical grid along the 2-surface u> = 0. This is reflected in 
the properties of the solution of the LYE. For the physical region we know that there is a unique bounded solution 
smooth up to ,y for which we cannot specify any data. In contrast, in the outside region we have to specify boundary 
data on the grid boundary in order to obtain a unique solution. This solution will also be smooth up to uj = but the 
fact that we can specify arbitrary data suggests that the two solutions will in general not match smoothly across J 1 . 
The physical reason behind this phenomenon is simple. We have two asymptotically flat space-times which are forced 
to match 'at infinity' in such a way that uj = is the J? + for one space-time and for the other. By suitably 
rescaling the metrics we can achieve that the two J^'s are isomorphic but there is no reason why the radiation which 
leaves one space-time should be the same as the radiation entering the other. Or, put differently, the mass-aspect 
obtained from one side is not necessarily the same as the one obtained from the other side. This 'mismatch' manifests 
itself in a jump in the third derivative of <j). While we cannot influence the part of the solution inside J 1 we can use 
the boundary data on the grid boundary to change the character of the space-time beyond . Thus, with general 
boundary data on the grid boundary we obtain a global solution which agrees with the unique solution inside J^, 
which is smooth for w/0 but which has a jump in the third derivative at uj = 0. 

This discontinuity will be present also in the initial data computed from the conformal factor Q and will propagate 
along J 1 like a shock front moving through a fluid. This is of course not what we are aiming for. What we need are 
initial data which are smooth on the entire grid. To achieve this we need to find a 'preferred extension' beyond J? of 
the unique solution from the inside. This is determined from the inside by the only requirement of smoothness. In the 



1 unless one would work with more than one grid 
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case of analytic data this is the unique analytic extension. The example of the Schwarzschild resp. Kruskal space-time 
given in the contribution by Schmidt fl3|| shows that in general we have to expect that there may be singularities of <p 
beyond but, hopefully, these will not be immediately adjacent to null- infinity so that in those cases which are of 
interest to us i.e., where the unphysical region is suitably small we can hope that the function (j) will be regular on 
the entire grid. 

So what can we do to obtain a sufficiently smooth solution of the LYE across J^? Unfortunately, this is not so easy 
and until now not solved in a satisfactory way. The first task is to obtain an idea about the boundary data which are 
necessary for the smooth extension. Fortunately, this is straightforward. The degeneracy of (Q) yields information 
about the behaviour of the solution on Setting to — in (Q) gives 

n a n a = -<j>\ (10) 

where '=' means 'equality on and where we have defined n a = d a uj. Furthermore, assuming that (f> is smooth we 
take a derivative of (0) and then put u — 0. This yields 

-2 (n e n e 5 a c + \n a n c ) d a <t> + (d a n c - \h ac d e n e ) n a = 0. (11) 

Since the map 

x a >-> (n e n e Sc + in a ri c ) x a 

can be inverted this equation and (^) allow us to determine the solution cf> up to its first derivatives on .y . We can 
use this information to estimate the value of </> on the grid boundary. We arrange lu in such a way that the outer region 
can be foliated by the leaves lu — const, and introduce a vector field s a which is transversal to them. This could be the 
vector field normal to this foliation but this is not necessary at this stage. We extend two arbitrary coordinates (£,77) 
defined on J* to the outer region up to the boundary by requiring that they be constant along s a . In this coordinate 
system the grid boundary is given by an equation of the form lu = cj(£, rf). Now we have s a d a <f> — d<p/ du and we can 
estimate the boundary data ips by 

V ),t, rf) « 0(0, f , r,) + |^(0, t rf) ■ lo(C, rf). (12) 

Clearly, this estimate can be useful only in those cases where J? is close to the grid boundary. Furthermore, it depends 
very much on the arbitrary vector field s a which is used to set up the isomorphism between J? and the grid boundary. 
There is an obvious choice for it: the vector field which is normal to <# with respect to the freely specified trial metric. 
Another possibility is to choose a vector field which is singled out by the numerical setup and which is transverse to 
both ^ and the grid boundary. 

With these estimated boundary data we can now proceed to solve the LYE on the entire grid. Numerically this does 
not pose a problem. For demonstration purposes we implemented a simple finite difference procedure to solve this 
boundary value problem using Octave ^l). Assuming axisymmetry we have the following ansatz for the free metric 

h ab = -f {dx 2 + dy 2 ) - g 2 r 2 dtf, (13) 

where / and g are arbitrary functions of (x,y) and r 2 = x 2 + y 2 . It is easy to convince oneself that the general 
axisymmetric hypersurface orthogonal metric in three dimensions contains two free functions. Therefore, the conformal 
class of these metrics is characterized by only one free function. Strictly speaking the ansatz above is too general. We 
find it rather useful to include this redundancy because it allows us to test the code on exact solutions. In Fig. |l| is 
shown one example of a solution of the LYE where we have specified </> = 1 on the boundary. There is nothing special 
about these boundary data except that they are positive. The solution is obtained on the square [0, 1] x [0, 1]. The 
functions /, g are assumed even in both x and y. The boundary defining function is u> = r 2 — (4/5) 2 so that J 5 " is at 
r = 0.8. The diagram on the left shows the solution itself, while the diagram on the right shows the third differences 
in the y-direction to illustrate the jump obtained in the third derivative of the solution. 

In Fig. H is shown the solution of the LYE with the same free data except for the boundary data which have been 
obtained from the estimate above. It is clearly visible that the jump in the third derivative has been greatly reduced 
(note the different scales in the two plots) but is still present. As an illustration we also show in Fig. [|the difference 
of the two solutions. Here it is clearly visible how effectively the LYE 'shields' the inside from the outside. The 
difference of the two solutions in the inside is several orders of magnitude less than in the outside. 

What can be done in order to find a smooth extension of the solution beyond tf! One procedure which comes to 
mind is the following: We start the solution process by estimating the boundary data and solve the LYE with these 
data. Once we are in convergence we switch from this 'inner' iteration to a different, 'outer' iteration. This time we 
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FIG. 2: Solution (left) of LYE and its third y differences (right) for estimated boundary values 



discretize the equation up to the boundary in the sense that we not only impose the equation at inner points of the 
grid but also on the boundary by using one-sided differences. Then we do not need any boundary values. Instead of 
prescribing the boundary data in this iteration we fix the value of the solution at some point in the interior which we 
know from the previous iteration. Fixing an interior point is necessary in order to prohibit the solution to converge 
to zero which is a valid solution of the LYE. Once this 'outer' iteration has converged we switch back to the 'inner' 
iteration prescribing the boundary data we have obtained from the outer iteration. Repeating these steps several 
times one hopes to obtain an increasingly smooth solution. This method has not been fully tested and it remains to 
be seen to what extent it can yield useful results. 

Another possibility to proceed is to use spectral methods. The reason for this is the fact that these methods are 
global on the grid and there is some hope that this globality helps to obtain a globally smooth solution. What follows 
is very much work in progress and so far there are no results in this direction. However, it is worthwhile to discuss 
the possible merits of this approach. 

Spectral methods are described in detail in various monographs || [ll], [H| Q so we can restrict ourselves to the 
essentials. The general idea is to use two representations of the grid variables / simultaneously. The first one is the 
grid representation where / is given by its values // = /(£z) at some grid points £j while the second one is obtained 
by writing / as a finite sum /(£) = Yli=o fl^liO where the functions Pj(£) belong to some set of 'basis functions'. 
The collections {fi}, the coordinate representation, and {fi}, the generalized Fourier representation can both be used 
to represent the function /. The transformation from one set to the other can be performed by multiplication with 
the matrix M l m = Pi(£ m ) or its inverse. The basis functions are chosen according to the specifics of the problem like 
symmetry, boundary conditions etc. For some specific functions like the trigonometric and Chebyshev polynomials 
the transformation from one representation to the other can be speeded up by applying fast transforms like the FFT 
in various forms. 

In the present problem we look for solutions of the LYE which are even on [—1,1] X [—1,1], restricted to the first 
quadrant [0, 1] x [0, 1]. Thus we may conveniently choose the tensor basis {Ui m (x, y) := T2i{x)T2m{y) '■ < l,m < N} 
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for some N. The functions T m are the Chebyshev polynomials of degree m defined on [—1, 1] by the relation 

T m (cost) — cosmt, for t 6 [0, 7r]. (14) 

The LYE can be written formally as 

Lcf) = ~6(j) 5 (15) 
where L is a linear differential operator with coefficients which vary with the coordinates (x,y). Inserting the ansatz 

N 

4>{x,y) = ^2 (j>imUi m (x,y) (16) 

l,m=0 

into ([is]) we find that we need to evaluate the derivatives of the Ui m and multiplicative terms of the form 
f(x, y)Ui m (x 1 y) for some functions f(x,y). Here, the advantage of having both representations at hand is appar- 
ent. We can evaluate the multiplicative terms in the coordinate representation while the derivatives are efficiently 
and accurately evaluated in Fourier space because there exist simple recurrence relations between the coefficients of 
a function and those of its derivatives (see any of the references cited above). Having evaluated these terms and 
transformed the result to Fourier space we can regard the linear part of the Yamabe operator as a linear operator in 
Fourier space represented by a 27V x 2/V-matrix L mapping the coefficients 4>i m to the coefficients of Lcf>. 

The nonlinear equation is solved by Richardson iteration (see e.g. (37)). We compute an update S<fi from a given 
estimate (f>r n \ of the solution by solving the linear equation 

L 5$ + 304>f n) 64> = - [U [n) + 6$\ n) ) (17) 

for 5(f) and use <f>r n +i) = 4>{n) +S(f> as an improved estimate. However, as it stands the matrix L is in general degenerate. 
When solving the equation using finite difference methods the same problem appears: the matrix of the discretized 
equation is singular. The remedy is, of course, to add boundary conditions to select a unique element in the kernel. 
But it is exactly these boundary conditions which are not known in our problem and the prescription of which leads 
to the discontinuity of the third derivative. 

The same applies in Fourier space. In order to get a unique solution we need to prescribe additional conditions. One 
aspect of these additional conditions must be their non-locality in the coordinate representation i.e., we should not 
try to fix any particular values of the solution at grid points. Such conditions would translate into global conditions in 
Fourier space. This means that we should impose local conditions in Fourier space e.g. by fixing certain coefficients. 
This affects the solution as a whole. 

A problem which arises is that in order to prescribe conditions in a reasonable way we should at least know how 
many. So far this is not quite clear. It is quite straightforward to convince oneself that the defect of (the spectral 
equivalent of) the Laplace operator obtained in the present setting is N + 1 and not 2N + 1 as one would have 
thought on first sight. The reason is that in the expansion only even polynomials occur and so the operator maps 
even polynomials to even polynomials, effectively halving the number of degrees of freedom. However, L contains not 
the Laplace operator but that operator multiplied by lo together with various other terms and at least the vanishing 
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of lj could enlarge the kernel. Numerical experiments show that this might be the case but they are inconclusive so 
far and a rigorous analysis is needed. 

Even if we figure out how many and which conditions to impose we need to be careful in applying them. For, 
suppose we solve the linearized equation subject to some linear condition B(8<t>) = then the nonlinear solution 
will necessarily satisfy B(cf>) — B{<p^). Thus, the condition will be 'remembered' in the nonlinear solution and this 
implies that we have to impose conditions for the linearized equation without imposing (restrictive) conditions for the 
nonlinear equation. Candidates for such conditions are those which eliminate higher degree polynomials. The reason 
is that one expects an exponential decrease in size of the spectral coefficients of the solution so that these should not 
contribute too much. 

An entirely different approach is based on a higher order approximation. Recall that the Richardson iteration 
scheme is obtained from the following consideration. Let N be a nonlinear operator (between finite dimensional 
spaces for our purposes) and assume we want to solve the (nonlinear) equation N[<j>] — 0. In order to derive the 
iteration procedure we obtain the equation for an update from an earlier estimate <f> by putting N[<p + Sip] = and 
linearizing with respect to the (presumably small) update 8<p. We get 

N[<j)]+dN[(t)] -84> = (18) 

from which the Richardson scheme follows by solving for Sip. Here, dN[<j>] is the linearization of N around (p. Suppose 
that dN[<j)] has a non-trivial finite kernel spanned by {Jia)a=i.... ,n for some n. Now we have the problem discussed 
above that we do not know what to do with the additional freedom in choosing the homogeneous solutions. However, 
we may take the expansion one step further 

N[<j) + 8<j)} w N[4>] + dN[4>] ■ S(f> + ^d 2 N[p] {8(f), 8<f>) . (19) 

Suppose Sep is a solution of (|is|). Then so is Sep + c a^a- Inserting this into jl9| ) we obtain 

N[<f> + 80 + ^ c A h A ] = \d 2 N[<P] (S</>, 8cf>) + £ c A d 2 N[<P] {8<f>, h A ) 

! A (20) 

+ 9 Yl c ACBd 2 N[(j)} (h A , M ■ 
1 A,B 

This equation can be considered as a quadratic equation for the coefficients ca and we cannot in general expect that 



this will have a solution. However, we can expect to find a set ca so that the right hand side of (20) is minimal. 
The procedure then would be to 

1. take an estimate (p of the solution 

2. compute the linearization dN[(j)) of the nonlinear operator TV 

3. determine a particular solution of flis} ) 

4. determine the kernel of dN[<fr] 

5. compute the second variation d 2 N[4>] and find the set of ca which minimizes ( |20| ) 

6. finally put 0^0 + c a^a and repeat from step 1. until convergence 

Steps 3. and 4. can be performed by finding the singular value decomposition p9| of dN[(j)]. Step 5. is easy in the 
case where N is the discretized LYE because then the second variation of N is simply 

(8!<p, 8 2 <f>) ' * 12O0 3 (5x0) (8 2 <t>). (21) 

This procedure provides a unique update and therefore a unique solution without imposing any conditions 'by hand'. 
Therefore, one could hope that this will lead to a sufficiently smooth extension of a solution of the LYE onto the 
exterior region. However, numerics is full of surprises and it is never known whether a theoretically devised procedure 
performs satisfactorily in practice. Work on this approach has just begun and it remains to be seen to what extent it 
is useful. 
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IV. CONSTRUCTING INITIAL DATA 



Now we turn our attention to the next problem which arises in the construction of initial data for the conformal field 
equations. As we discussed above (cf. sect. II) we have to obtain the initial data for the Ricci and Weyl tensor once the 
correct conformal factor which selects the metric with constant negative curvature from the given conformal class has 
been found. The problem we have to face is the fact that we need to divide two quantities which vanish somewhere. 
It is guaranteed from the analytical results that these quantities vanish simultaneously so that the quotient is well 
defined and (at least) continuous. However, numerically at least one of the quantities is affected by numerical error 
which implies that it will not vanish at those places where the other quantity has its zeros. 

We write any one of the components in (JsJ) or (||) as q = C/lo where C contains the specific component and the 
factor (j) 2 which comes from the conformal factor fl in the denominator. Then we know the denominator (the freely 
specified boundary describing function ui) exactly and its computation at one grid point involves only round-off error. 
On the other hand, the numerator C can be computed only up to a certain accuracy which determined essentially by 
the truncation error of the discretization of the LYE. So we have C — C + e. We denote by q = C/oj resp. q = C/w 
the computed and exact results of the division of the computed resp. exact numerator C resp. C. Let G = {xfc} 
denote the grid and its points on which we perform the computation. Then we have 



\q(xk) - q{x k )\ 



Hxk)\ 
K^fc)|' 



(22) 



The maximal deviation from the exact result occurs at the grid point where |w(a;fe)| attains its minimum. This can be 
arbitrarily small, even zero, if u> happens to vanish at a grid point. If |e(xfc)| would be only due to round-off error as 
it would be if we computed the exact analytic solution then this division would not be disastrous (if we take care to 
prevent the vanishing of ui on any grid point). This is because division is a backward stable process in IEEE floating 
point arithmetic [^9| [fl], [f9|. Then we divide two non-zero numbers which are 'simultaneously' small. The result is a 
well defined number which approximates the correct quotient within round-off accuracy. 

The real disaster occurs when e is due to the truncation error. Then e(x) may be orders of magnitude different from 
the minimum value of u> so that the quotient will be wrong by orders of magnitude. In this situation the numerical 
implementation of l'Hopital's rule to compute the indeterminate term 0/0 in one form or another is of no use. The 
reason is that this rule essentially replaces the differential quotient by a difference quotient. This provides a good 
approximation if one uses the exact expressions for the numerator. However, if the numerator had been previously 
computed with a truncation error then we will also get an answer which is off by orders of magnitude. To illustrate 
the problem let us consider a simple example (cf. Fig. [|). Take u> = 4/5 — x 2 on [— 1, 1] and let q be some smooth 
function and define C — qto. Thus, C is known exactly and can be computed up to round-off errors. The diagram on 
the left shows the error \C/u> — q\. Clearly, this is on the level of machine precision. The peaks are not related to the 
location of the zeros of w indicated by the small triangles on the x-axis. 




-0.5 0.0 



FIG. 4: Division by u with (right) and without (left) additional random perturbation of the numerator 



Next we 'spoil' C by a random distribution of numbers with maximal size below 10~ 10 to mimic the effect of the 
truncation error (although truncation error is certainly not random) and compute again the error. The result is shown 
in the diagram on the right. Here, we see that the overall error is on the level of the random perturbation. At the 
location of the zeros the error is maximal. Thus, we see from this example that it is the structure of l/u> on the grid 
which determines the peaks. 

There have been two ways trying to avoid the disaster. Hubner |33| ] proposes to solve a singularly elliptic equation 
for the quotient in terms of numerator and denominator. The equation has the same degeneracy on ^ as the LYE 
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and hence one has to expect the same problems. In |M| we considered a procedure based on spectral ideas to perform 
this division. We want to elaborate on this method here. The main idea, summarized schematically in Fig. a is as 



C(x) 



FCT 



Ca 



•U(X 



IFCT 



FIG. 5: The main idea involved in division by u> 



follows. We represent the numerator and the denominator both as a finite sum of basis functions, say Chebyshev 
polynomials (the following discussion applies just as well to Fourier series). Multiplication by ui is a linear map which 
we can represent in Fourier space as a matrix lj 1 j defined by 

uj(x)-T j (x)='£uj i j T i (x). (23) 

i 

Given the representation of C in terms of its coefficients C^ then one could simply solve the linear equation q l = of ' j& 
to obtain the coefficients of the quotient. However, there are some snags. 

The basis functions {Tj : < j < N} span a /miie-dimensional vector space Vn- Since multiplication with any 
polynomial (except To) increases the degree we should obtain products with degree higher than N. However, if we 
regard the multiplication map as a map from Vn into itself we are asking to represent a high degree (> N) polynomial 
as a sum of polynomials of degree up to N. With this 'folding back' of high degree polynomials into the range of low 
degree polynomials - this effect is known as aliasing - multiplication is no longer injective. Therefore, we need to 
enlarge the target space. Multiplying two polynomials of degree at most iV yields a product of degree at most 2N. 
So we regard multiplication by as a map uj : Vn — > V%n- This map is injective so that the image w[Vat] of Vn is a 
N + 1-dimensional subspace in Vzn- The matrix u> of the map uj is no longer square but has dimensions (2iV + 1) x N . 

Let & be the representation of C — J2j=o in Vjv which is trivially also a representation in Vzn- In general, & 
will not lie in w[V)v] because the error in its computation will drive it away. Hence, the first step in the computation 
of the quotient is to project C onto w[Vat] and then the second step is the inversion of w on its image. 

Both these steps can be accomplished by computing the (reduced) QR-factorization of uj (cf. p9[). This algorithm 
allows us to write the matrix uj uniquely as uj = QR with a square (N + 1) x (N + 1) upper triangular matrix R and 
a (2N + 1) x (N + 1) matrix Q whose columns are orthonormal. The crucial fact for us is that these columns span 
w[Vat]. This implies that P = QQ l is a projector onto that space while Q t Q = id,N is the identity on Vjv- 

Thus, to solve the equation q = ujC we proceed in two steps: first compute z = Q t C and then solve q = 
These steps correspond exactly to the two steps mentioned above because Qz = PC is the projection of C onto w[Vjv]. 
Since R is the matrix representation of uj with respect to the basis provided by the columns of Q and its pre-image 
in Vn its inverse provides that vector whose image under uj is the projection of C onto a; [Vat]. 

The scalar product involved in the computation of the QR-factorisation is the natural scalar product between 
the Chebyshev polynomials. Hence the projector is the orthogonal projector onto w[Vat] with respect to this scalar 
product. The projection PC is that vector in w[Vjv] which is closest to C with respect to that scalar product. This 
is not what is needed here. For suppose that w is a polynomial of degree M < N. Then the projection should not 
contain any polynomials up to degree M — 1 because any of these would give rise to a singular quotient. Thus, we 
need to construct a different projection P which has the same image as P but which annihilates all polynomials with 
degree less than M. 

In order to do this we have to impose the condition that w be a polynomial with degree M < N. This is not a 
severe restriction because uj can be specified freely. Then we consider uj ; Vn — > Vn+m as a map into Vn+m- To 
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find the new projection it is necessary now to compute the full QR-factorisation of oj because this provides not only 
information about the image of uj but also about its orthogonal complement. In fact, one can again write uj = QR 
where now Q is an orthogonal (N + M + 1) x (N + M + 1) matrix and R is upper triangular with dimensions 
(N + M + 1) x (N + 1). Again the first N + 1 columns of Q form an orthonormal basis of w[V/v] while the remaining 
columns span its orthogonal complement, the kernel K of P = QQ* , so K = kerP. We denote by Qi the matrix 
containing the first N + 1 columns of Q and by (§2 the remaining ones. Then Q — (Qi\Q2)- Note, that Qi coincides 
with the matrix Q obtained above from the reduced QR-factorisation. 

Let Vm be the space spanned by the polynomials with degree less than the degree of lo. Clearly, we have 

V M nuj[V N }={0}. (24) 

We need to change the projector P = QQ l — Q\Q\ into a new projector P = Q\Q\ + Q1SQ2 for some (N + M) x M 
matrix S which has to be found from the fact that P annihilates Vm- This form of the new projector follows from 
the fact that in the adapted basis given by the columns of Q a general projector with image wfVjv] has the matrix 
representation 



liV 







, 



Let U be a (N + M) x M-matrix whose columns span Vm- Then S has to satisfy the condition 

(Q\ + SQl) = 0. (25) 



Hence we get 



Note, that Q 2 U is a square matrix of dimension M x M which is necessarily invertible. This is a consequence of (|2 



S = Q\U(Q\U\ 1 (26) 

and consequently 

(Qluy'Ql)- (27) 



P = P I 1 N+M - U 

Using this 'improved' projection kills all the lower degree polynomials and hence yields much smoother quotients. As 
an example consider Fig. ^| where we have illustrated a one-dimensional case with and without the use of the improved 
projection. The upper diagrams show the results obtained by dividing the function / = ui(2x 3 — 7x 5 + 10cos 2 (5a;)) 
by lo = (x — 1/2) (x + 2/5) (x — 1/100). However, to simulate the real problem we have contaminated / with a low 
order polynomial (i.e. at most quadratic in this case) with random coefficients of size at most e = 10~ 2 . The left 
diagram shows the quotient (indicated by a line with crosses) without the use of the regularizing projection. The line 
is the exact result. Clearly, the quotient is influenced by the random deviation. The difference to the exact quotient 
clearly shows the structure of 1/lu. The maximum error is influenced by two sources: the size of e i.e., the size of 
the truncation error in the application we have in mind. The second influence is the resolution, i.e. the number of 
grid points (or, equivalently, the number of polynomials). With higher resolution 1/lo can be better resolved which 
results in higher spikes and hence in bigger errors. On the right hand side the regularizing projection has been used. 
Here, the deviation from the exact result is not visible and, in fact, the error is at the level of 10~ 12 . Here, the error 
is not affected by e because the projection exactly kills any low order polynomial. It appears that the regularizing 
projection redistributes the error in a global manner so as to obtain a function which is 'evenly' divisible by to. 

As a final remark it should be noted that the way we have obtained the regularizing projection does not in any 
way depend on the dimension of the problem although we have essentially treated only the one-dimensional case. 
In a similar way one can treat higher dimensional cases by employing tensor bases and interpreting the spaces of 
coefficients appropriately. The geometry behind the scenes remains exactly the same. 



V. CONCLUSION 



As we have seen the question of getting smooth initial data beyond the conformal boundary of space-time is difficult. 
At the moment the only way to obtain these data is via the Andersson-Chrusciel-Friedrich procedure which has its 
limitations. Therefore, we should look for other ways to get such initial data. There are two other possibilities for 
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this. One way is to understand the structure of the conformal constraints and to find ways of formulating well-posed 
boundary value problems for their direct solution. This is the route taken by Butscher |j| . 

Another possible approach towards the construction of smooth initial data is suggested by recent work by J. 
Corvino |Q (refined by Chrusciel and Delay (T2J ) who has constructed smooth solutions to the constraint equations on 
an asymptotically Euclidean space- like hypersurface which coincide exactly with Schwarzschild data outside a compact 
set. The evolution of such data will contain at least for some time an asymptotic region which is exactly Schwarzschild. 
Therefore, it is possible, at least in principle, to find hyperboloidal hypersurfaces in the time development of the data 
on which hyperboloidal data are implied which are exactly Schwarzschild outside a compact set and these data can 
trivially be extended analytically beyond J? . To obtain these data it is necessary to either generalize Corvino's 
method to the hyperboloidal setting and solve appropriate equations directly there or to generate the hyperboloidal 
data numerically from Corvino's asymptotically flat data. This requires an appropriate evolution scheme for Friedrich's 
regular finite initial value problem at spatial infinity pa which so far is lacking. 

Another issue which ultimately has to be thought about is the inclusion of matter fields into the conformal frame- 
work. Including matter fields is complicated because the Bianchi equation for the rescaled Weyl tensor acquires a 
source term and this involves derivatives of the stress-energy-tensor. So far there has been only one case in which 
a system of matter fields coupled to the conformal field equations has been studied numerically |30|]. While there 
is no problem (at least in principle apart from complexity) to include fundamental conformally invariant fields like 
the Maxwell, Yang-Mills fields [^2| or other massless fields like the conformally coupled scalar field |n| problems can 
arise for massive fields like the Klein-Gordon field or ideal fluids. The reason is that the matter equations are not 
conformally invariant so that a regularization of the equations on ^ by rescaling the fields appropriately is unlikely. 
However, in the usual scenario of an isolated system in which a source generates waves and sends them out to infinity 
the matter region will not extend out to infinity so that the singular behaviour on is not present. Even if one 
assumes that matter fields extend to infinity one can still hope to be safe on J? because in some cases like the massive 
Klein-Gordon field the time evolution forces the field to decay exponentially towards J? so that the formal singularity 
is in fact not there |5(J. 

To conclude it is probably safe to say that the conformal approach has reached the status where we know that it is 
a reasonable and feasible approach which offers some exciting possibilities but which also provides us with problems 
(mostly of a numerical nature) which are (apart from some idiosyncracies) very similar to those encountered in other 
approaches. A joint effort of the numerical relativity community seems to be called for in order to overcome these 
problems. 
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